LARGE DEFORMATION ANALYSIS OF 
THERMO-ELASTIC CONTACT PROBLEMS 


by 

Utsa Majumder 



Department of Aerospace Engineering 
Indian Institute of Technology, Kanpur 


June, 2005 



Th 

i^t\%oos Im 

1 1 StP 2005 

*,#=rrT ^ 

VotW 

«nf<l (. A..!^'-'- 



Certificate 


This is to certify that the work contained in the thesis entitled Large Deformation Analysis of Thermo- 
elastic Contact Problems submitted by Utsa Majumder is carried out under my supervision. 

This work has not been submitted elsewhere for a degree. 




Dr. C.S. Upadhyay 

Department of Aerospace Engineering 
I.I.T. Kanpur 


Dedicated to: 


My Parents and Teachers 



Acknowledgements 


First of all, I would like to convey my gratitude to my thesis supervisor Dr. C.S. Upadhyay, whose 
able, professional and yet friendly guidance and support was with me for more than one year as a 
source of motivation. His interests and confidence in me were the underlying source of inspiration for 
this work. 

Moreover, 1 am thankful to all my friends who made my stay at IIT Kanpur a memorable one. Spe- 
cially, my days with Hall-7, E top members will be a sweet memory for the rest of my life. Some of 
them, like Shamik-da, Krishnendu, Anik, Soumya, Partha, were very friendly and helpful through out 
my program. I would also like to name some of my friends who are not in the institute now, and some 
seniors like Ashok, Pritam, Somnath, Abhijit-da, Sidhhartha-da, Sandip-da, Kaustubh-da, Rakesh-da, 
Kaushik-da, Debasis-da and many others who shared their happy moments with me. 

1 shall like to thank Abhijit-da specially for providing some helpful technical hints related to my 
work. Also the help, support and motivation from every person related to this work directly or indirectly, 
are gratefully acknowledged. 


Utsa Majumder 
June, 2005 



Abstract 


The present work is mainly targeted to develop an algorithm for total and updated Lagrangian formu- 
lations tor large detormation dynamic thermoelastic problem with contact from their basic definitions 
using proper transformations and assumptions. 

The main trame-work of this study is based on the Updated Lagrangian formulation and its application 
for some standard problems by coding it using Finite Element approximations. Updated Lagrangian 
method is developed on the incremental basis from the current configuration to the unknown deformed 
one. Related objective stress measures are formulated using the required incremental transformations. 
With proper approximations, the Jaumann-Zaremba objective stress increment rate is used for the com- 
puter program. 

Thermal equations are considered to be decoupled from their elastic counterparts and are solved inde- 
pendently beforehand to generate required temperature fields and thermal stresses which is a input to 
the elastic problem. 

Simple contact condition like sticking contact with a rigid wall is considered. The contact constraints 
are resolved into discrete form required for Finite Element Analysis using the Lagrangian multiplier 
method. 

Also some work had been done on continuation methods and advanced solution procedure like the 
arc-length methods. The existing theories are revisited and a specific method, namely the cylindrical 
arc-length method, is incorporated into the static code to capture a few well known unstable phenom- 
ena. 
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Chapter 1 
Introduction 


The primary goal of the study is to understand the actual meaning of the still debatable terms - updated 
and total Lagrangian methods. An idea, clear to some extent, has been developed from the existing 
materials. The basic definitions of the processes are used to develop the corresponding balance laws. 
The resulting equilibrium equations are expressed in an incremental form. This use of incremental 
form introduces the requirement of the use of objective stress measures to confirm frame indifference 
of the physical parameters. The existing objective stress measures like Jaumann-Zaremba stress rate 
and Tmesdell rates are shown as the special fomi of the objective stress increments and are used with 
proper approximations. 

In this study, decoupled thermo-elastic problem is solved. By this approximation, first the temperature 
field in the body is obtained by solving the heat conduction problems,using finite elasticity formulation. 
This temperature field is then used to obtain the themial loading on the structure. In the scope of this 
work, though the constitutive theories are deduced in a general manner, the balance equations for heat 
flux are considered to be decoupled from the elastic part while developing the computer code. 

The contact problem, one of the most attractive and yet difficult one, is considered for a simple sticking 
contact with a rigid body. Basically, the challenging problem of multi-body deformable mechanics 
with contact search problem is reduced to a mere single body problem with pre-imposed constraints on 
certain degree of freedoms. 

Also some continuation techniques are studied regarding the advanced solution methods. These limit 
point analyses are very good examples where the updated method has an positive advantage over the 
total Lagrangian method. Instead of handling the intricacies of the complex continuation techniques, 
rather a simple method like cylindrical arc-length method which is a readily programmable one, is 
introduced in the code. 

Ultimately, for the dynamic analysis of the models, Newmark p integration scheme is used for the 
time integration part. Though this is an implicit time integration scheme, it shows some oscillations 
when implemented to the code. This phenomenon is also mentioned in some of the literatures. The 
oscillations resulting from the standard integration scheme proposed by Newmark has been damped 
out by introducing algorithmic damping. 


1 



1 . 1 Literature review 


As the topic of the thesis was a classical one, many text books are available regarding the issues. 
Though some of them give conflicting definitions on updated Lagrangian method, overall they all focus 
on the same procedure and tactics. 

Gurtin [1] and Malvern [2] has given a nice treatise on classical continuum mechanics which is essen- 
tial for the understanding of the elementary kinematics of the deformable bodies. The basics of updated 
Lagrangian scheme are stated in Bathe [3] and Belytschko [4]. Both of them mentioned the updated 
and total Lagrangian approach as the extension of the same kind of analysis. Belytschko [4] has de- 
scribed the updated Lagrangian method to be the one where the measures of the physical quantities and 
the domain of the integral are taken about the Eulerian coordinate instead of the Lagrangian one used 
in total Lagrangian formulation and also showed the equivalence of the two formulations by proper 
coordinate transformation. But Bathe [3] has marked the difference between these two by the updation 
of the initial configuration used. According to him, in updated Lagrangian method, the initial config- 
uration is updated continuously by one step to be used in the next step. The advantage of the updated 
Lagrangian method, as described here, lies in the predicting the actual path for large deformation prob- 
lems (application of Taylor’s series expansion is more valid as the total deformation path is divided into 
small parts) as well as the numerical advantages. Though Bathe, Ramm and Wilson [5] mentioned the 
incremental equations of motion, the expression of the conservation laws were in total or instantaneous 
forrn and hence, for the kind of transformation used, no objective increment in the stress tensor was 
done. Gadala and Wang [6] deduced an incremental approach for analyzing the kinematic behavior of 
deformable bodies which they described as the ALE (Arbitrary Lagrangian Eulerian) approach. This 
approach closely follows the incremental formulation used in this thesis. The concept of the intermedi- 
ate reference frame used in the ALE formulation by Askes, Kuhl and Steinmann [7] is also very close 
to the intermediate reference frame at time t used in this analysis with the original frame at time 0 and 
current frame at time t + At. McMeeking and Rice [8] presented an Eulerian approach to tackle large 
elasto-plastic deformation problem which, when transformed to the reference configuration, produces 
the Lagrangian method. 

Crisfield [9] has given some applications of the updated Lagrangian method. He also introduced some 
objective measures of stress increment rate [10] while treating large strains. The objectivity require- 
ments for different kinematical quantities are mentioned in [1 ],[1 1 ] and [2]. Lecture notes on mechanics 
of solids of Brown University available from the web [12] have elaborately deduced the Oldroyd stress 
rate as a general objective rate for large deformation and large rotation and has shown that the other 
stress rates are special versions of the Oldroyd rate. Belytschko [4] also mentioned about the Jaumann, 
Truesdell and Green-Naghdi rates according to the constitutive equation used. With relevant assump- 
tions, the Jaumann-Zaremba stress rate is used in the computer program as shown in [10]. 

The thermal conduction laws and governing equations for heat transfer are available in Bathe [3], Cook 
[13] and Haupt. Haupt [11] and Lubliner [14] also treated the constitutive assumptions for thermo- 
elasticity in a very general way. The constitutive assumptions mentioned in [1] and [2] are used to 
develop the constitutive relations for the analysis. In the code, isothermal material constants, given in 
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Lubliner [14], are used with standard numerical values. 

The contact formulations used in this study are devoid of any kind of complexity. Rigorous study on 
the contact conditions for both small and large deformation is done by Wriggers [15] for sticking and 
slipping contact. Several schemes like Lagrangian multiplier method, penalty method, perturbed La- 
grangian method etc. to handle the constraints are discussed thoroughly. The constraint conditions are 
discretized to fit the finite element analysis following the Lagrangian multiplier method and consider- 
ing the contact to be a sticking one. Bathe [3], Crisfiield [10] and Belytschko[4] have also given some 
brief review of the contact mechanics. 

In this study, some work has been carried out on the advanced solution algorithms, mainly focusing on 
the continuation techniques. This methods generally guides the solution after a limit point is reached. 
A short textbook review of those methods are included in Crisfield [9, 10]. Similar works and their 
applications have been carried out by Ahmed and Xiao-Zu [16]. The cylindrical arc-length method 
introduced by Crisfield [9], is incorporated in the code and are applied to some standard structures with 
limit point loading. 

For dynamic analysis of structures, implicit time integration scheme has been used. The schemes, 
explicit and implicit are available in Bathe [3] and Belytschko [4]. Newmark-/? scheme, proposed 
originally by Newmark, is used in the computer program. It is simple to understand, use and apply in 
coding. But, the problem with this scheme is that, though it is said to be unconditionally stable, very 
often it shows an oscillation about the equilibrium values. Cook [13] mentioned the unstable or slow 
convergence behavior of Newmark integration scheme and has given some conditions on the time step 
length for the same to be unconditionally stable. This issue is also addressed by Mahato [17]. Some part 
of his work was regarding the stability of the Newmark-/? scheme and he also introduced algorithmic 
damping by increasing the values of the Newmark constants. He mentioned the term arrested insta- 
bility about this issue and also addressed slow convergence, though in [4] it is mentioned with respect 
to the explicit time integration scheme. In the current implementation, some oscillations are observed 
in the solution obtained from the program and these oscillations are damped out increasing the values 
of the Newmark constants to some higher numerical values. Further study using a more sophisticated 
time integration scheme, should be carried out in order to resolve this issue. 

1.2 Structure of the thesis 

In the current study, the work and simulations carried out have been addressed in an arranged man- 
ner. Chapter 2 covers the description of the different reference systems used in continuum mechanics, 
transformation between them, definitions of Eulerian, updated and total Lagrangian methods, basic 
conservation laws stemming out of the Reynold’s transport theorem etc. The theory of objectivity has 
also been addressed in this section. This chapter is a short overview of the tools used in the whole study 
and are presented from a general point of view. The basic weak forms are developed and incremental 
techniques are described here. 

In the next chapter, the focus is on the formulations done in this study. The continuum mechanics 
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laws and transformations presented in Chapter 2 are applied to the deformable solid body model to 
get the incremental version of the weak form of momentum balance. Using the laws of objectivity, 
the Jaumann-Zaremba rate, which is used in the current formulation while discretizing the equilibrium 
equations, is shown as a special case of the general objective stress rate with some assumptions. 

The constitutive assumptions are also mentioned in this chapter and the constitutive laws are developed 
initially for hyper-elastic materials and then for thermoelastic materials using related potentials. 

In Chapter 4, the contact constraints are formulated and discretized, Lagrangian multiplier technique is 
used and most of the formulations are done following Wriggers [15]. 

Chapter 5 contains overview of some continuation methods, namely the arc-length methods. Spherical, 
cylindrical and linearized arc length control techniques have been discussed under the light of general- 
ized load-displacement control. 

Chapter 6 is about discretizing the weak forms of conservation equations obtained in the previous chap- 
ters. Different strain displacement relations are formulated for Lagrangian shape functions. Newmark- 
0 time integration scheme is used to modify the differential form into an algebraic equation. 

Chapter 7 contains the results obtained from some simulations done using the computer program and 
their comparison with existing commercial codes and/or strength of material formulation. Mainly dif- 
ferent static analyses are emphasized in the chapter. The issue of numerical oscillation for dynamic 
problems are also addressed briefly. 
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Chapter 2 

Preliminaries of Continuum Mechanics 


2.1 Kinematics 

The target of this chapter is to focus on the tools used to analyze the systems under consideration. 
For kinematic systems, the basic laws are governed by continuum mechanics. The laws of nature here 
means the balance laws that the body must obey, namely, the continuity equation, conservation of linear 
and angular momentum and conservation of energy. In continuum mechanics formulations, generally 
two coordinate systems are used to completely describe the kinematics of a body as shown in Figure 
2. 1 for the convenience of the analysis [ 1 ]. The space fixed (with respect to the stars) coordinate system 



Figure 2.1: Reference, undeformed and deformed coordinate systems 


is referred as the Eulerian or the spatial coordinate system and the one moving fixed with the body (or 
sometimes fixed with respect to some intermediate position) is called the Lagrangian or the material 
one. The governing laws of mechanics are written with respect to any one of these two reference 
systems. And correspondingly the type of the analysis changes. There are mainly two types of analysis 
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done in mechanics, Eulerian or spatial formulation and Lagrangian or material formulation. These two 
different formulations are just the different representation of the same balance laws with respect to the 
two different coordinate systems mentioned above. 

In Eulerian or spatial formulation, the equilibrium equations are written with respect to the spatial 
coordinate system. The discretization is also done in the Eulerian coordinate. And for Eulerian finite 
element formulation, spatially fixed meshes are used. In contrast with this, in Lagrangian formulation, 
the equations are referred with respect to a material fixed coordinate. The basic idea of Lagrangian 
formulation is to fi'eeze the reference coordinate system in time. For static problems, it is frozen to 
some virtual time step. Updated and total Lagrangian methods only differ in the regard of fixing the 
coordinate to a different time step [3]. 


2.2 Transformation between the two coordinate systems 

The domain of the body in Euclidean space, Bt is property by which the bodies are referred from both 
material and spatial configuration. As shown in Figure 2.2, let us assume that the material coordinate 
of the body occupying the domain Bt is expressed by X and the corresponding spatial coordinate for it 
is X. 



From the basics of continuum mechanics, the spatial description of the body is a mapping / of the 
material coordinate (and vice-versa) as; 


5 = Kx) 


( 2 . 2 , 1 ) 


The expression detVf represents the local point-wise volume after deformation. And / is assumed to 
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be a one-to-one mapping of the points X, it implies that the body cannot penetrate itself [1]. So, 


detVf > 0 

The vector u{X) is the displacement between two generic points in the spatial and the material refer- 
ences. 


f{X) =X + u{X) 

The tensor quantity, defined point-wise, 

^ = V/(X) 
dx 

~ M 

is called the deformation gradient tensor [1]. As detVf > 0, deformation gradient F is a positive 
definite tensor. Deformation gradient may be 
dx = FdX 

F = 


Mathematically, the deformation gradient is the Jacobian matrix of the motion of the body x{X, t) [4]. 
It includes both rigid rotation and stretch. The total deformation may be considered as a rigid rotation 
operating on a stretch or vise-versa. Hence, polar decomposition may be used to find out the stretch 
part separately from the rigid rotation, which causes no strain with respect to the material reference 
ffame[l]. 

F = RU (Right polar decomposition) 

= (Left polar decomposition) 

where, R is the rigid rotation tensor, U and V are the right and left stretch tensor respectively. 

As R is an orthogonal tensor,the Right Cauchy-Green Strain tensor and the Left Cauchy-Green Strain 
tensoA^ defined as: 


: considered as a measure of deformation or strain as 


dx 

M 

dJiX) 

dX 

-?.(X + u) 
dX^ 

- du 
1 + —^ 

8X 
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^ ^2 zT z 

Right Cauchy-Green strain tensor; C = U = F F 

z z2 Z zT 

Right Cauchy-Green strain tensor: B = V = FF 


Hence, the expressions for right and left Cauchy-Green strain can be derived as [1] 

: zT z 

C = F F 

dX dX 

~ du ,duj , du J , du . 

= l + -^ + i-^) +{—^) — 
dX dX dX dX 

And, B = pf 

dX^^ dx’ 

~ du ,duj , du du P 

The tensors U and V are similar, where C and B are also similar as; 

P = Mr 

z z z zT 

And, B = RCR 


where ^ is an orthogonal matrix such that R ^ = R'^. 


2.3 Balance Laws 

2.3.1 Reynold’s transportation theorem 

The balance laws of continuum mechanics, or the conservation laws, are all branching out of the trans- 
portation theorem [11, 1], Transportation theorem is basically a balance relation for an open system. 
The general statement of the Reynold’s transportation theorem is as follows: 

Let $ be a smooth spatial field, and assume that $ is either vector valued or scaler valued. Then for 
any part Bt and time t. 


d_ 

dt 

£ 

dt 


I 

JBt 

I 

JBt 






4> -I- ^diVx V ) dQ 


/.( 

/ / 

■iBt Jd. 


■ ndS 


(2.3.1) 
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Proof of Equation 2.3.1 is available in Gurtin [1] and Haupt [1 1]. This famous equation is basically 
the origin of all the balance laws of nature. A general open system is shown in Figure 2.3. The region 



Figure 2.3; Control volume for time window At 

7Z is the region under consideration for control volume analysis. For control mass analysis, which is 
followed in this section to get the corresponding conservation laws, TZ is equivalent to Bt at time t and 
to Bt+At at time t + At. 


2.3.2 Continuity equation 

The most important property of the bodies those are defined in the Eulerian space, is the matter con- 
tained in that space of the body B or in other words, mass. Mass of a body is defined as: 


m(X) 


/. 

L 


dm 

ptdQ 


Using Reynold’s transport theorem, 

d 

—m 
dt 

The term ptv ■ MS is a surface term representing the out flux of mass through the surface. But as 
Bt defines the body as a whole, this term is zero. Hence, 



/ p'tdQ. -f / ptV- ndS 

JBt JdBt 
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Now, if the deformation can be represented as a mapping / and its gradient F = V/ then, transforming 
domain of integration of L.H.S to Fo [1], 

det (#(A0) Pt{x{X))dn = j poiX)dn 

Hence, the continuity (of mass) equation turns out to be: 

det (#(A)) ptix) = po{X) (2.3.2) 

Let $ be a continuous spatial scaler field. Then, for any given body Bt, 


j ^{x)pt{x)dD. = 

JBt 


^X)det{F)pt{X)dn 


I Bo 



^x)po{x)dn 


So, as a point-wise statement, it can be written that. 




^X)po(x)da 


(2.3.3) 


2.3.3 Momentum balance 

The mechanical motions of a body is governed generally by the external causes named ’’Force”. In this 
study, the forces may be categorized mainly as: 

(i) Contact forces between different parts of the body or the internal forces, 

(ii) Contact forces exerted on the boundary of the body by the environment or external surface forces, 
and 

(iii) Forces exerted on the volume of the body by the environment or external body forces. 

The most important axiom in mechanics to describe the balance of momentum is Cauchy's hypothesis. 
This axiom assumes that there exists a surface force density function T(x,t) for each unit normal 
vector h and every (x, t) on the boundary dBt as shown in Figure 2.4 [1]. 

This contact force is often referenced as the traction force. For the whole boundary dBt, the total 
traction force is 

/ f{n)dS 

JdBt 


The sum of the body force applied by the environment on the body can be expressed as 



boundary 


Figure 2.4: Surface force density function 


Hence, the statement of the balance of linear momentum becomes: 

Total external force, /eit(X,f) = j f{n)dS+ j 

JdBt JBt 

By Newton’s second law, if the linear momentum of the body is £, then 

dC 


bdfl 


(2.3.4) 


fext — 


dt 


(2.3.5) 


And by Reynold’s transport theorem, 

d£ d 


dt 


pvdfl 


= -f 

JBt 

= f ^(pv)dQ,+ f {pv)v ■ hdS 
JBt dt JdBt 

= / pvdQ, + / vpdQ, 

JBt Jst 


'Bt 

as u • n is a measure of flux through the boundary dBt and is zero for control mass formulation. 
Hence, 

dC 


dt 


— mu + mu = fext 


(2.3.6) 


From Equation 2.3.4 and Equation 2.3.6, as the mass of the system is time independent, the force 
balance equations turns out to.be 


/ f(n)dS + f bdQ, = j pudQ, (= mu) 

JdBt Jst Jst 


(2.3.7) 


Cauchy’s theorem (existence of stress): [1] 

This theorem defines a spatial tensor quantity a, called the Cauchy stress tensor. The existence of the 
stress tensor is a necessary and sufficient condition that the momentum balance equation should satisfy. 
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The statement of the theorem is; 

(i) For each unit normal n, 3 a tensor <5- such that, f{h) = dh 

(ii) a is symmetric, and 

(iii) a satisfies the equation of motion 

divia + b = pu (2.3.8) 

The last statement is a point-wise statement of momentum balance and can be derived by putting the 
expression of a from statement (i) in Equation 2.3.7 and using Gauss divergence theorem. 


2.3.4 Theorem of virtual work 


Theorem of virtual work (or theorem of virtual power for the corresponding rate fomi) is the direct 
outcome of the momentum balance. A kinematically admissible (obeying the essential boundary con- 
dition) virtual displacement function Su is assumed from the space all possible displacement fields of 
the body. The point-wise momentum balance equation (Equation 2.3.8) is multiplied and integrated 
over the whole domain to get the balance of virtual power. 


/ [divxd' 4- b] • 6udQ. = pu - 8udCi 

JBt JBt 

or, j divi [a) ■ SudQ. ■+■ pb- 6udQ. = u - 6udCt 

JBt JBt JBt 

using integration by parts, 

or, / [divx [aSu) — a : gradx5u]dQ 3- / pb - Sudfl — j pu - Sudil 
JBt JBt JBt 

using Gauss divergence theorem, 

/ dSu - ndS —I a : qradxSudO. -t- / pb - 8ud€l = / pu - 6ud^ 
JdBt JBt JBt JBt 

using statement (i) of Cauchy’s theorem. 


/ 

JdBt 


T - 6udS 


/ a : gradxSudO, + j pb ■ 6udQ, = j pu - 6udQ, 
JBt -JBt JBt 


So, the theorem of virtual work or the weak form of the momentum balance equation in the deformed 
coordinate is: 

f pu ■ 6udD, 3- /" a gradxSudO, — j pb - 8udQ, — T - 5udS = 0 (2.3.9) 

JBt JBt -JBt JdBt 


2.3.5 Balance of Heat Flux 

The governing heat conduction equation is build over the assumption that the material imder consider- 
ation follows the Fourier’s Law of heat conduction. Then for the heat flux equilibrium for any internal 
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point of the body [13, 3] 


dxi 


‘-Xl 


de 

dxi 


0x2 


36 


A 

dxz 


36 


^^3xz 


+ q 


de „ 


(2.3.10) 


Assuming an admissible temperature field 56 , the theorem of virtual work done by the heat flux, can be 
derived as: 



=>■ [ {ki6^i)^idQ, + f g^SddQ — [ cp^56dQ = 0 

Jst JBt JBt et 


+ - cp—\59da = 0 

dt 


Hence, the weak form turns out to be 



[Sd{ki9^i)]nidS 


[ {Ki9,i)5e^d^+ f 

JBt Jb 


iH9dQ, 


Bt 


L 


cp969dVL = 0 


(2.3.11) 


Balance of heat flux and corresponding weak forms are more elaborately discussed in Section 3.5. 


2,4 Change in observer and Invariance of reference frame 

Use of two coordinate systems, like the spatial and the material ones, or the co-rotational one, which 
moves fixed to some specific direction with the body, is one of the artistry of mechanics that eases the 
analysis of the system. Use of two reference coordinates, or rather two different observers generally 
means that there exists two different views of the body with respect to the different observers. Let us 
consider two motions i(t) and x*{t) of a body Bt as shown in Figure 2.5. Then x and x* is related by 



Figure 2.5: Two reference systems 
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the relation 


x*{X, t) = q{t) + Q{t)[x{X, t) - 0] (2.4.1) 

for every material point X and time t, where g is a point in space and Q is a rotation. This is equivalent 
to define a one-to-one mapping /(i, t) as; 

J{x,t) = g(0 + - 0) 


and representing x* as the deformation x followed by a rigid deformation /' . That is, 

As shown in Figure 2.5, the two observer coordinates x and x* are related as: 

X* = f{x, t) 


And the inverse mapping. 


~ r — 1 ~ SIC 

X = J ox 


is also tme. Inverting this relation to the material coordinate frame, it turns out to be 


X{.,t) = PU)of{.,t) 


The term objectivity means the invariance of the various kinematic quantities under the change of 
reference system. The laws of mechanics should strictly follow the transformation laws between the 
two deformation mappings x and x*. The basic relations of the different quantities are derived in this 
section. 

Let, 

dX 

- dx” 
and F* = -^ 
dX 


differentiating Equation 2.4.1, 

F*{X, t) = Q{t)F{X, t) (2.4.2) 


* As both the frames x and x* are mutually orthogonal triad, any one of them can be transformed to the other simply by 
a rigid rotation, namely Q 
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Equation 2.4.2 gives the transformation law for the deformation gradient. As the tensor Q represents a 
rigid deformation, det{Q) = 1, i.e. 

det{F) = det(F’) [as det(ab)=det(a) det(b)] (2.4.3) 

This confirms the tact that the volume of the deformed body Bt should remain same for both the 
observers. 

Also, let the polar decomposition of the deformation gradients be 

F* = R*U* = V*k 

I' = Ml = ^h 

From Equation 2.4.2, 

F* = ku* =QP 
= Qru 

Since both Q and R are elements of the vector space formed by the rotation tensors, 

k = Qh , (2.4.4) 

and 

U*=U (2.4.5) 

i.e the stretch U remains same for the two different observers. The only difference in the deformation 
gradient comes from the rotational parts. The corresponding rotation should be transformed as in Equa- 
tion 2.4.4. From Equation 2.4.5, the right Cauchy-Green strain tensors for the two different observers 
also remains same 

and, d* = 

= u*^k^kd* 

= U*'^U* [as R* is orthogonal] 
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And the left Cauchy-Green strain tensors are related as follows: 


B* = F*F*^ 

= [as = P] 

Now, referring to Figure 2.6, the axiom of change of observer states that, 



Figure 2.6: Surface force density function with respect to the two different observers 


ft* = Qn 

And as 

We can also write 

f*{h*) = Of (n) 

Using T = an and T* = a*n*, 

f* = d*n* = Qf 
= Qan 

= QaQ'^h* [as 0"’^ = Q^] 

Hence, 

d-*{x*,t) = Q{t)a{x, pCP^it) (2.4.6) 

Equation 2.4.6 gives the transformation rule for the stress tensor for a rigid rotation Q to the observer’s 
reference system. 
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2.5 Eulerian approach 


As mentioned earlier, the basic idea of the Eulerian type of analysis is to study the continuum with the 
observer fixed at the Eulerian or spatial coordinate system. It is like studying the kinematic and kinetic 
characteristics of some spatially fixed (with respect to the stars) points for some interval of time. As this 
approach, by definition, accepts that the material points under study (at the fixed point at current time) 
may not, or rather, shall not be the unique one, so, it is not possible to study the motion or displacement 
of any material point in this scheme. Instead, the velocity and accelerations are considered to be the 
primary parameters of investigation at the fixed points of interest. This understanding is somewhat 
analogous to the control volume analysis of fluids. 

2.6 Lagrangian approach 

Here, the idea is that the parameters under study are the ones related to the material points. The 
formulation is in terms of Lagrangian measures of stress and strain, and the derivatives and integrals 
are defined over the material reference frame. The material coordinates are expressed as some discrete 
function of the spatial coordinates. This is obtained by freezing the time at the instant of obseiwation. 
The instant at which the system is frozen to the observer defines the types of Lagrangian methods, 
namely the Updated Lagrangian or Total Lagrangian. 

In Total Lagrangian(TL) method, the configuration of the body at time t=0 (fi-om when the analysis 
starts) is taken as the reference configuration for the whole analysis. The balance laws of mechanics 
are transformed to the configuration at time 0 from the currently deformed, unknown configuration. 
Whereas in Updated Lagrangian(UL) approach, the reference configuration is continuously updated 
with time, i.e., the whole time period under analysis is divided into some small time windows and 
for each time window, the TL approach is carried out. For the first step, the reference configuration 
is the body at time 0. And for the next step onwards, the reference configuration used in a specific 
UL step is nothing but the deformed configuration (output firom the previous step) obtained from the 
previous step. Both UL and TL formulations are theoretically equivalent, but the only difference for 
UL method is that the configuration at time t is considered as the initial configuration, where for TL, it 
is the configuration at time 0 [3]. Hence the UL method is nothing but the TL method by steps and is 
used while handling non-linear problems to get the numerics better. 

2.7 UL Method and incremental form 

The target in UL method, as aforementioned, is to solve for the unknown parameters over some pre- 
defined time steps, called UL steps, each of time size At {At may vary fi'om step to step). For any 
generic time step with initial time as t, the unknown parameters are to be solved for the time t + At 
from the known configuration at current time t. As the configuration of the body at that time, Bt+At is 
an unknown itself, it is not possible to solve the momentum balance equation directly for time t + At 
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(as it is basically a force balance over the unknown domain). The individual terms of the momentum 
balance law as well as the domain of definition of it, must be transformed to the ones at currently known 
time by using proper transformation laws. For small deformations, these two configurations, and 
Bt may be considered to be the same and the the momentum balance equation or its weak form can be 
solved directly in a total form defined over the domain at time t. But for large deformation and large 
rotation systems, this may lead to a completely different solution than the exact one. 

But the problem in this transformation, basically transformation of the equation of the principle of vir- 
tual work from the deformed configuration to the undeformed one, makes the resulting equation a rather 
complex one and with a lot more nonjinear terms. And in this transformation, it is also convenient to 
use the incremental form instead of the total form. Detailed deductions are provided in Section 3.2. 
This is nothing but defining the primary unknowns at time t -I- At as a linear increment over their current 
values as: 


■tit+At = Ut + Ati (2.7.1) 

and 

o’t+At = dt + Act (2.7.2) 


and so on. 



For the static analysis, the quasi-static form of the UL method is used. As there is no time parameter in 
the formulation for static problem, the total force Fext is divided in some incremental forces like AFi, 
AF 2 and so on as shown in the Figure 2.7. The system is analyzed for each force step cumulatively for 
each AFext. The configuration obtained from the previous is the input configuration for the next step. 
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Chapter 3 

Theoretical formulations 


3.1 Balance law in deformed configuration 


In the current work, the deduction of the momentum conservation equation and corresponding weak 
forms are done by giving proper honor to the instances of time when they are defined and hence to 
the related reference frames also. The governing differential equation(GDE) is the momentum balance 
equation defined at the deformed configuration. The conservation of momentum, defined at a specific 
configuration, say Bt, is expressed as: 

[ fdS+ f pbdn= [ pudQ. (3.1.1) 

JdBt JBt JBt 


Replacing T from Cauchy’s theorem asT — ah, and using Gauss divergence theorem, the momentum 
balmce equation can be written in a point- wise sense as: 


/ [diva ph](Kl = / pudil 
JBt JBt 


(3.1.2) 


Perturbing the system by a kinematically admissible virtual displacement field 6u the weak form of the 
momentum balance becomes analogous to the theorem of virtual work (Equation 2.3.9) in the deformed 
configuration Bp. 


f pu-6udn+ f a:.^dn- f pl-Suda- [ f-6udS = 0 (3.1.3) 

JBt JBt ' 9X Jst JdBt 

3.2 Transformation of weak form to undeformed coordinate and 
incremental weak form 

The target of Lagrangian formulation is to bring the momentum balance equation from the unknown 
deformed configuration to the currently known- and undeformed configuration. As shown in Figure 
2.1, the momentum balance equation is defined at the configuration Bt+M which is currently unknown. 
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So, the target is to form the momentum conservation law, or rather it’s weak form in the deformed 
configuration and then transfer the same to the undeformed configuration for that time step (i.e to the 
Lagrangian reference frame defined over time t). This transformation from t+Atto time t is dependent 
on the displacement gradient from time t -t- At to t, which is shown as tensor G in Figure 2.1. 

This displacement gradient G, may be called the incremental displacement gradient as it is analogous 
to the other incremental parameters from time t to t -I- At, and can be found out in terms of unknown 
incremental displacements from Figure 2.1, using the definition of the corresponding parameters as 
follows. 


and 



dxt 

dX 


ax 


d{X -f Ut+At) 
d{X + fit + Afi) 


■■Ft + 


dX 

dAu 


So, 


Ft+At — { I + 


dAu 

dxt 


(3.2.1) 


And by definition of the incremental displacement gradient, it operates on the displacement gradient at 
time t(.Ff)to give the displacement gradient at time t + At(Ft+At)- That means, 


h+At = doP, 


(3.2.2) 


Comparing Equation 3.2.1 and Equation 3.2.2, the incremental deformation gradient G can be ex- 
pressed as: 



(3.2.3) 


To transform the weak form from the deformed to the undeformed coordinate, first of all, the weak 
forms at those configurations are to be found. And for that, the kinematically admissible displacement 
fields at both time instances t and t + At should be known. But as the essential boundary conditions 
are also defined at time t, hence, only the admissible displacement field at time t is known. 

Here, it is assumed that for both time instances, the displacement field has the same kinematic con- 
straints. So, it is possible to use the same virtual displacement Su for both the cases. 

Then, the momentum equilibrium at the deformed configuration Bt+At is just the same like Equation 
3.1.3 with the integration limit changed over the domain of the body at time t + At (i.e the deformed 
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configuration), 


pu • SudVL + 




f z dSu 

/ O' : -T^dO, 

ox 


-i 


pb ■ 6udQ> —I T • 6udT = 0 


(3.2.4) 


Transforming Equation 3.2.4 to the configuration at time t and writing it in point-wise form [1], 

dSu 


iB, 


d&tG\div : jdfi -|- / Pt^t+At ■ dudVt — / Pt'^t+At ■ dudG (3.2.5) 

0^ JBt Jst 


And the weak form of momentum balance at time t itself is, 

f dSu f f 

/ [div (dtSu) — dt : -1- / pibt ■ SudQ, = / ptUt ■ 5udQ (3.2.6) 

Jbi Ox 


Eliminating Equation 3.2.6 from Equation 3.2.5, the incremental form of the virtual work turns out to 
be, 


I [detG {div - dt+At ■ ~ {^tSu) + dt : 


+ / Pt[bt+At - it] • SudO. = / pt[fit+At - ut] ■ 
JBt JBt 


(3.2.7) 


6udQ, 


Assuming Ut+At = fit + Afi and fit+At = fit + Afi, Equation 3.2.7 can be modified to have Afi as the 
primary unknown in it as. 


/ [detG div [dt+At^u) - div [dtSu)]dO, + / ptAb-Sudfl 

JBt JBt 


j. 

JBt 


dSu 

PtAu • 5udQ, 4- / Ad : -^dfl 


L 


(3.2.8) 


dx 


The first term of Equation 3.2.8 can be expanded as, 

/ detG div {dt+At^u) — div (dtSu) dQ,= detG dt+AtSu ■ ht — dtbu • 
JBt ^ JdBt 


nt 


IdBt 


detG dt+Apk ■ bu — dtfit ■ Su 


dS 

dS 


IdBt 


detG {dt + Ad)nt ■ Su — dtUt ■ Su 


dS 


Hence, 


detG div {dt+AtSu) — div {dtSu) 


dCl 


IdBt 


detG {dt 4- Ad)nt ■ Su — dtUt ■ Su 


dS (3.2.9) 


It is assumed that the traction forces on boundary do not change with configuration of the body, i.e 
there are no follower loads. Then, the effect of change of normal on the traction force term shown in 
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Equation 3.2.9 can be thrown out. 

Thus, on applying the Gauss divergence theorem, 

/ ptAi ■ 5udn +[ Ad: = f ptAh ■ Suda + I AT- 6udS (3.2.10) 

JB, JBt ox 

Equation 3 .2. 1 0 is the incremental weak form of the momentum equation in material coordinates. 

The assumption that the kinematic admissible conditions on the virtual displacement field 5u remains 
same for both configurations Bt and Bt+At, is valid as there is no change in the kinematic constraints in 
the two configurations. 

Obviously, this assumption doesn’t remain valid for large deformation contact problems. Then the 
constraints are very likely to be changed within the time span At continuously. For those kind of for- 
mulation, the virtual displacement field should also be considered as a function of the spatial coordinate 
and time, i.e (Ju = t). 

3.3 Objective stress increment - Jaumann Zaremba stress rate for 
small deformation 

As mentioned in Section 2.4, objectivity is an issue of importance when two types of reference frames, 
namely spatial and material reference frames, are considered in the analysis. The kinematics of a 
body appears to be different to two different observers situated in the two different reference flames. 
The basic idea of using the objectivite measiures is nothing but just honoring the difference between the 
observations of the two observers. Objectivity for rigid rotation generally means a set of transformation 
mle to be implemented between two reference coordinates to transform some quantities defined in one 
coordinate to another. 

Equation 2.4. 1 and Equation 2.4.6 state the general rules for objective transformation of any kinematic 
vector and kinematic tensor respectively for the configurations shown in Figure 2.5. 

For large strain and large rotation analysis, the value of the increment in Cauchy stress due to the rigid 
rotation of the deformed coordinate with respect to the undeformed one, is never negligible. For the 
incremental displacement field Au, the rigid rotation part of the whole deformation from time t to t+At 
can be found out by decomposing the corresponding incremental displacement gradient G = {I + 
by polar decomposition. Let the decomposition of the incremental displacement gradient be as, 

d = M 

by the right polar decomposition. 

Here, R is the rigid rotation tensor and U is the stretch tensor. 
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As R is an orthogonal tensor, 


= G^G 

dAu / dAu\^ /dAu\''dAu 

^ dxt ^ \ dxt ) \ dxt J dxt 


Assuming the incremental deformations considered to be small and using binomial expansion 


u = yu^ 


/ ~ dAu f dAu 

~ y V 

~ If dAu dAxJ 

~ / + 1 

2 \ dxt dxt 

= !+b 


where D is the symmetric part of the tensor and 


k = d(j^ 


dxt J \ 2 V dxt dxt J J 


= ^ 1 / dAu dAuF 


2 V dxt 

= i + W 


dxt 


where W is the rigid rotational part of the tensor 

The tensors D and W are understood to be some measure of linear and rotational strains respectively. 
Hence, it can be stated that the transformation G can be equivalently decomposed in two cumulative 
deformation gradients namely, the stretch U and rotation Ft. The increment in stress will also have two 
corresponding increments. The pure stretch part D causes the stress increment as 

i^Ktr.tcH='ED ( 3 . 3 . 1 ) 

^ As any vector can be considered as the Eigen vector of the isotropic tensor /, hence Eigen directions of D and 4- , 

D being any general tensor, are same and the second one can be expanded using Binomial theorem. 
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where E is the corresponding material tensor, and the rigid rotation R causes the stress tensor in the 
deformed coordinate to be 


rotation RcTfR 

= {! + W)at{i + 

= dt + Wat + a-WatW'^ 

The higher order term WatW^ is neglected’ to linearize the objective stress increment Aarotaticm as 

i^~^)rotation = ^'^t + at^^ (3.3.2) 

From Equation 3.3.1 and Equation 3.3.2 the total stress increment turns out to be 

Ad = Id + (3.3.3) 

By definition, the stress increment rate a can be written as 


d’f+At — + ED + W^t + 



a = ED + Wat + ^tW^ 


The Jaumann-Zaremba stress rate is defined as the objective increment rate of Cauchy stress and is 
given by: 

= ED = a - (4^at -h atW^) (3.3.4) 

This stress rate is the objective rate of Cauchy stress that comes out of the theoretical formulations 
' automatically when the analysis is carried out in the spatial reference frame and assumption of small 
deformation is considered. The most general expression for objective stress rate, the Oldroyd stress 
rate [14, 12], and its derivation is mientioned in Section 3.4. Jaumaim-Zaremba stress rate is a derivative 
of Oldroyd stress rate that takes care of only the linear terms. 


3.4 More general objective stress rate and related approximation 

Refer to Figure 2.1, the incremental deformation gradient G is decomposed into a rigid rotation and a 
stretch. Expressions for these two components are derived in Section 3.3. 

'The discussion on such throwing out of higher order terms and keeping the quantities still objective at the same time, 
is given in Section 3.4 
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For large strain and / or large rotations, maintaining objectivity for stress increments is very much 
important and this enforces the transformation mles already mentioned. In Section 3.3, some higher 
order terms are intentionally thrown out of the formulations to linearize the equations. In objective 
formulations, neglecting some higher order term may invite considerably large errors unless the terms 
thrown out .is objective itself. Before neglecting any term in this kind of calculations, one must be 
careful so as not to destroy objective requirements. 

3.4.1 Linear expression for objective stress increment 

For a rigid rotation R, the objective stress transformation is given by 

d* = Raff 

and as shown in Section 3.3, this can be expanded as 

a* = a + wa + aw^ + waw^ 

From the definition of objectivity, the second order term WatW'^ is itself an objective quantity as W is 
the orthogonal tensor which represents the rigid rotation between the two coordinate systems. So, the 
first order frame indifference approximation can be made by dropping the second order term, which 
keeps the rest objective also as Aa = W at atW'^. This gives exactly the Jaumann-Zaremba stress 
increment rate for Cauchy stress without destroying the objectivity axioms up to the first order terms. 

3.4.2 Oldroyd stress rate 

As shown in Section 3.3, the Jaumann-Zaremba stress rate inherits cancellation of non-linear terms. 
But the derivation of the most general objective stress rate is quite different. The Oldroyd rate, some- 
times referred as the general objective rate [12], is rather easy to be deduced by the rate formulation. 
The transformation fi'om one reference frame to the next ’star marked’ reference fi:ame with deforma- 
tion gradient as G is considered. Cauchy’s stress is a second order tensor and transformed as 

a* = PaP'‘^ (3.4.1) 


Differentiating both sides. 


= PaF^ + PaP'^ + PaP'^ 


Using Equation 3.4.1, 


I* = 


a* + PaP'^ + r P-^P'^ 
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As the velocity gradient L can be expressed as L = FF ^ [ 1 , 11], the obj ective rate of d* can be written 
as 

d* = + F (3.4.2) 

Equation 3 .4.2 gives the material frame invariant relation for the time derivative of Cauchy stress tensor. 
The Oldroyd objective stress rate is defined as 

Om. = 5 — L a — a U" (3.4.3) 

This objective rate is sometimes considered as the general objective stress rate [12] and the Jaumann 
stress rate is deduced from it considering rotation followed by stretch (or the reverse), i.e L = W. 
But the process of polar decomposition seems more general as the two motions are separated as two 
operators. 

3.5 Heat transfer 

In the current analysis, the coupling between the effects of elastic and the thermal loading is ignored. 
That is, the thermo-mechanical constitutive equations (Equation 3.6. 13) are considered to be isentropic. 
This is basically ignoring the second law of thermodynamics and removing the constraints on the en- 
tropy generation. When the analysis is assumed to be isentropic, having small quasi-static incremental 
steps, this kind of formulation gives a simpler system to solve. 

The material under analysis is assumed to obey Fourier’s law of heat conduction. Then the heat flux 
in the body in the direction is given by [3] 

qi = - kij9j {i,j =lto3) 

For a three-dimensional body, the general form of the balance of heat flux is [13] 

de 

ikxx9 kxyO^y -j- kxz9 ,z\x (^kyx9^x kyyO^y -f" ky^O d” ikzx^ ,X d” kxyO^y F kxz9 ^z) Q ^ 

(3.5.1) 

As the medium is considered to be isotropic and homogeneous, = 0 and Equation 3.5. 1 becomes 

(^^X^.x).X d- {kye,y),y F (L^..,),. d* - CpO = 0 (3.5.2) 
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Multiplying Equation 3.5.2 by a thermally admissible virtual temperature field SO and integrating over 
the domain of the body, 




{^i0t-\-At,i)4S6dil^ + 




qlAt^oan 


cpt-^/:^t0t+AtS0dQ = 0 








/ i^i^t+At^i^S) ,id^l — I (^k.i9t^^if)S6^idri 

I ^Pt-\-At0t-{-AtS0d^ 0 

+ At 

Using Gauss divergence theorem, the weak form at deformed configuration becomes 

I (,ki6f,^^i^-iS9')n-idS j (^ki6i^/\i.i)S9j,dQ,-h f q^_^_^^69dVl — f cpi^^t9t+Atd9dQ: = 0 

(3.5.3) 

Weak form of Heat flow equation in the previous configuration is 

/ {ki9i^i89)nidS - f {ki9t,i)d9,idi'l + [ q\89di\ - f cpt9tS9dQ = 0 (3.5.4) 

As in Section 3.2, the increment in temperature is also considered a linear one as 9t+At = + A9. 

Eliminating Equation 3.5.4 from Equation 3.5.3, assuming the higher order terms from direction change 
of surface normal (as shown in Section 3 .2) to be negligible small and ignoring the coupling with elastic 
deformation, the incremental weak form for heat conduction is 

I cptAe89dn+ [ {kiA9^i)S9^ida = f Aq'’S9dn+ f Aqlni89dS (3.5.5) 

JBt JBt ' ' JBt JdBt 


'Bt JBt JBt 

where qf is the surface heat flux in the direction. 


3.6 Constitutive law and corresponding assumptions 

The axioms of nature of forces, or rather the laws they should follow, are defined by the momentum 
conservation equations. These do not define the deformation completely. Because, the relation between 
the amount of force and the corresponding amount of deformation is still unknown. To relate the defor- 
mation and it’s type to the force field, one has to apply required constitutive assumptions to the body. 
In mechanics, there are mainly 3 types of constitutive assumptions used [1], namely, 

(i) Constraints on the possible deformation field of the body. This assumption defines the t 3 q)e of defor- 
mation that would take place for the force. Example of this kind of constraint is like allowing the rigid 
motions of the body or imposing an constraint on the volume of the body like isochoric deformation. 

(ii) Assumption in the form of the stress tensor. This kind of assumptions defines the properties of the 
stress tensor, e.g. symmetry for no body moments, diagonal in case of pressure load etc. 

(iii) Constitutive assumptions relating stress to deformation. These are the conventionally used consti- 
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tutive equations like Hooke’s law. They express the measures of force as a function of the motion. 

The constitutive assumptions are basically the statements of relation between the displacement field and 
force fields, consistent with the kinematic process C (C is called the constitutive class of the body[l]) 
which enforces the compatibility conditions. 

3.6.1 Theories of constitutive equations 

The constitutive relations must obey some laws evolved from the fundamental postulates of mechanics. 
Most important of those are [2, 11] 

(i) principle of equipresence and 

(ii) material frame indifference. 

Principle of equipresence 

Though this axiom has not been universally accepted, but this has a very clear physical sense and is 
followed by many of the analysts. This principle was first presented by Truesdell and Toupin (1960) 
as: 

An independent variable assumed to be present in one constitutive equation of a material should be 
assumed to be present in all constitutive equations of the same material, unless its presence contradicts 
the principle of material frame-indifference or some other fundamental principle [2], 

Material frame indifference 

This assumption states that every constitutive equation should be objective, or independent of the rota- 
tion of the material frame they are defined at (frame indifference included in Section 2.4). This is one 
of the fundamental postulates of a purely mechanical body. 

3.6.2 Isotropic linear elasticity 

The assumption of traditional elasticity is that the stress-strain relation is not dependent on any other 
parameters. In other words, the general measure of tension can be represented as a function of only the 
general deformation [14, 11]. i.e. 

f = f{i) 

where f is the general tension parameter and e is the any general deformation. 

And also, in general, the word elasticity is analogous to hyper-elasticity. This constitutive class in- 
cludes materials which have a strain energy function defined due to Green as[14] 


W = W(l) 

These materials are also called Green-elastic. Elasticity without any underlying strain energy function 
is referred as hypoelasticity or Cauchy-elasticity. 

Now, for small deviation in the deformation parameter, the constitutive relations may be linearized and 
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the strain-energy function W may be expanded about an equilibrium value Wq by Taylor’s series as 




^0 + 


de 


e+ higher order terms 


(3.6.1) 


Considering f to be the work conjugate of e, the work done about the equilibrium state is 

W(l) = Wo + rl 

Comparing Equation 3.6.1 and Equation 3.6.2, 


0 

£ = 0 


(3.6.2) 


0 di 




Hence, 


0 


T = ::r Ot time t=0 

de 


(3.6.3) 


Now, Generalized Hooke’s law states that each stress component depends on all of the strain compo- 
nents and vice- versa[ 14]. i.e 

Tij = GijkiSki (3.6.4) 

where Cyfc; is called the elastic constant. 

From Equation 3.6.3 and Equation 3.6.4, 


gw 

deij 


— ^ijkl^kl 


Or, 


Ciikl — 


d^w 


(3.6.5) 


— r\ o 

oeijOSki 

The constant Cijki is a 4th order tensor and is called the material tensor. There are 81 entries on that 
tensor, but from the symmetry 






deijdski dskidEij 

Cijki is symmetric, and only 36 terms are significant. 

There is also another additional symmetry for Cijki which can be shown by changing the indices as 



if the strain measure e is symmetric' Thus the number of significant constitutive constants become 2 1 . 
For isotropic materials, for directional symmetry, this number further reduces to 9. 

For the Cauchy stress and its work conjugate Almansi strain, the linear elastic constitutive equation is 

U’ij ' (3.6.6) 

Cij, the isentropic tensor of rank 4 can be expressed as [14] 


Cjjfcj — XSijSki + jJ'SikSji 4- k'SiiSjk 


Using the symmetry of Cijki and using Lame’s constants. 


— X5ij6kk T 


(3.6.7) 


Using engineering constants, namely Young’s modulus E and Poisson’s ratio i/, the relation of Equation 
3.6.7 becomes[14] 


a 


ij — 


E 

(1^ 


[(1 “F ^^kk^ij\ 


(3.6.8) 


For three-dimensional bodies, if the six independent entities of the Cauchy stress aij and Almansi strain 
€^is)written in the vector form following the Voigt notation [4] as dk and Cfc, then the material tensor 
becomes a two-dimensional tensor, which operates on e to give ct. 

As shown in [3], the material tensor then takes a form as^^^^ 


[C] 
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(I-U 

E{i - y) 

(1 + i/)(l - 2u) 0 

0 
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(i-U 
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(I-U 
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(1-U 
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(1-U 
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0 
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0 

l-2v 

2 ( 1 - 1 /) 

0 
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0 

0 

0 

0 

l-2u 

2 ( 1 - 1 /) 
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0 

u 

0 

0 

0 

l-2t/ 

2(1-UJ 


(3.6.9) 


3.6.3 Thermoelastic constitutive relations 


For thermodynamic constitutive law, r and e are the thermodynamic tension and the thermodynamic 
deformations respectively. 

The different potentials used for thermodynamics are defined as [1 1] 


The internal energy U is a potential for thermodynamic tension r for isentropic process and Helmholtz 
free energy is the potential for isothermal process. 

The relations among the potentials and the independent variable can be deduced from the thermo- 


'ey =e,,,hence^ = -fi^ 


dei 
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Thermodynamic potentials 

Relation 

Independent variables 

Internal Energy (Id) 

U 

5, Sij 

Helmholtz Free Energy (^) 

^ = U~s8 

^5 ^ij 

Enthalpy (h) 

h Id rijSij 

S,Sij 

Gibbs Function/Free Enthalpy (g) 

g = U - sd — TijEij = h — s9 

d,rij 


Table 3.1: Thermodynamics potentials 


mechanics as [ 1 1 ] 




(3.6.10) 


Using Cauchy stress cTy in place of and Almansi strain in place of £y , and expanding Helmholtz 

free energy by Taylor series about an equilibrium configuration 0[14], 




5^ 


'^>{6, €ij) = 4'o + L(^ ~ ^o) + ~ ^0 o) + higher order terms 


dtij 


(3.6.11) 


The relation between for the thermodynamic potential 'i’ can also be shown to be analogous to that of 
Equation 3.6.5 


^2^ 

deijOCki 


(3.6.12) 


Assumed the constant values 


de^ 0 

dOdsij 0 
^deijdcki 0 


- 7r<^o 


-/?° 


daij 


dcki 


= C 


0 

ijkl 


Now, from Equation 3.6.10 and Equation 3.6.1 1, the constitutive relations for linear thermo-elasticity 
can be shown as follows 


S = So + — ^o) + -Pij^ij 

VO P 

- ^ij ~ -60) + 
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For isentropic process, s = sq, hence Equation 3.6.13 reduces to 




POCQ 




= 4 + + -f 4/; 


pao 




(3.6.14) 


Using the engineering constants and the coefficient of thermal expansion a, the stress-strain relation 
stain becomes [14] 


tij = —[(1 4- u)aij - uakkSij] 4- a{e - 
Inverting Equation 3.6.15, the stress cry can be expressed as 


(3.6.15) 




(1 - 


(3.6.16) 


where C is the linear elastic constitutive coefficients and = (^ — ^o)- i-e, in Voigt matrix notation. 
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(1 4- i/)(l - 2z/) 
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2(1 -UJ 

1 0 0 0 0 0 ' 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 


aAd 


(3.6.17) 
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Chapter 4 

Analysis of contact constraints 


The contact problem is a very complex one among all the mechanical analyses presently. As computers 
are used to solve the huge systems of equations, study of contact has become more solvable and inter- 
esting. In this analysis, a two body contact is formulated. The hitting or slave body is assumed to hit 
a target or master body (Figure 4.1) [15], The general formulation for two deformable body systems 
is done in this chapter, though for the computational part, only contact with a rigid target body was 
considered. The main problem in analysis of large deformation contact analysis is that the assumption 


Initial configuration Deformed configuration 



Figure 4. 1 ; Initial and deformed configuration of contacting bodies - minimum distance 

made in Section 3.2 that the kinematically admissible displacement field Su remains same for both the 
configurations Bt and Bt+^t does not stand valid. Su should be considered as a fimction of time t and 
the spatial coordinate x also. 

But as far as small deformation contact problems are concerned, this assumption can be followed with- 
out a large amount of error caused to the results if the time steps are taken judiciously, i.e the time steps 
should be discretized in such a way that at the start of a particular time step, the bodies should be in 
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contact with each other. 


4.1 Normal contact of three-dimensional bodies 

In Figure 4. 1 , two deformable bodies are shown at their initial and deformed configuration. The contact 
analysis between two material points and X^, poses a minimum distance problem to be solved in 
the deformed configuration as [15] 


Pl -^2ll = -52|| (4.1.1) 

X 2 CI 2 

Having the point X 2 defined on the boundary r 2 from Equation 4.1.1, or in other words, dropping the 
normal on the master surface from the slave point (considering the master surface to be at least a locally 
convex region, as shown in Figure 4.2), the non-penetration condition for normal contact can be stated 
as 

{xi - X 2 ) ■ h'^ > 0 (4.1.2) 

For analysis of the contact problem. Equation 4. 1 .2 must be imposed on the standard weak form of the 
problem. 


4.2 Lagrangian multiplier method - penetration function 


Hitting body 



Figure 4.2: Initial and deformed configuration of contacting bodies - minimum distance 

Lagrangian multiplier method is the best way to get the actual solution from the constrained (Equa- 
tion 4.1.2) boundary value problem. For this approach, the non-penetration condition is modified in a 
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penetration function as shown in Figure 4.2 [15] 


g{xi,x2) 


(.ri - X 2 ) if (xi - i's) < 0 
0 otherwise 


(4.2.1) 


This functions is basically the gap between the two nearest points affer contact [15]. The are multiplied 
with the unknown Lagrangian multipliers and added to the potential energy of the system, when the 
corresponding Lagrangian multipliers A as the corresponding contact forces at those nodes. 


4.2.1 Sticking contact 

For this case, the work done by the contact forces are given by 


n 


cont 



(4.2.2) 


The contact work from Equation 4.2.2 is added to the kinemetic potential to get the total potential 
energy as 


ntotat ri/cm T ^cont 
or , "F ^^amt ~ 0 


This gives the weak form for the contact problem as 


/ pu ■ Sudi} + 
Jb 



- f pb- Sudn - ! f- 5udS - [(X-Sg + g- SX)dn = 0 (4.2.3) 

ox Jq JgB JB 


Again, using the assumption of small deformation that the kinematic constraints for admissible g and A 
remain the same for both the configurations Bt and Bt+At (i-e considering virtual increment in g and A 
in time t and t + At to be same as in Section 3.2), the weak form of momentum equation at time t + At 
can be written as 


[ Pt+At'dt+At ' dudO, + f dt+At ■ ^'d^ f Pt+Atbt+.At ■ budCl 




L 


Tt+At ■ SudS 

dSt+At 

/ (^t+At • ^9 + gt+At • SX)dQ, = 0 

d fit+At 

(4.2.4) 


Transforming Equation 4.2.4 to the undeformed configuration at time t and eliminating the weak form 
at time t itself from it as done in Section 3.2, the incremental weak form, after neglecting the higher 
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order terms, is 


PtAu ■ 6udQ +[ Ad: ^dQ, - f ptAb • 6udQ. - [ 
Jb, dX 


AT • 6udS 


(4.2.5) 


(AA • 6g 4- Ag.5X)di^l = 0 


4.2.2 Slipping contact with friction 


The formulation is same as the sticking one, only the contact potential energy can be expressed as 
follows instead of Equation 4.2.2 (as the constitutive assumptions in normal and tangential directions 
are different) 


{^nSn + • 9T)dS 


(4.2.6) 


where the subscripts N and T means the normal and the tangential components respectively. 

For frictional contact, Equation 4.2.6 is the constraint condition along with the corresponding frictional 
laws. If Coulomb’s linear friction law is used, the additional constraint is given by the constitutive 


assumption 



(4.2.7) 


where A^i and Aj-j are the two components of Ay. 

Now, let the tangent and normal at the point f j on the contact surface Tcont be as shown in Figure 4.3. 
The direction of the tangent is not unique, it can be taken along any convenient direction orthogonal to 


n’ 



Figure 4.3: Normal and tangent at the point of contact 


the direction of the normal. and ^2 are the curvilinear coordinate system at point X 2 on surface T 2 . 
Suffixes Ti and T 2 represents the components along tangent direction Ti and r 2 . 
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Let (f) be the angle between g and gr, . Then, from Coulomb’s friction laws, 

Xti = i^Xn cos ^ 

Xx2 — sin (j) 

Xti 

and, |A;^/| = -A^i seccf) 

And for discretization, it is also important to know the direction of the frictional force applied at the 
point of contact, By Coulomb’s friction law, it is in the reverse direction of the relative slip. 

Direction of relative slip between the elements at the contact is Aut and can be related to Au as 

Aux = Auxi + A'Ut’j 

= (Au ■ f,)fi + {Au ■ f2)f2 

= {f I ®fi)Au + (Tacgira) Au (4.2.9) 

= (Tj (Si Tj + 72 (8> T2 )Au 

= f'Au 

Equation 4.2.6, 4.2.7, 4.2.8 and 4.2.9 together gives the constraint conditions applied at the contact 
point for the frictional slip problem [15]. 

Then the rest of the formulation follows the same way as the no-slip one. 
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Chapter 5 


A brief study on continuation techniques 

Continuation techniques, as the name suggests, are aimed at tracing the complete path of the load 
displacement curve for a structure even after a point of singularity comes in the system. In general, 
when a limit point is reached, i.e the slope of the load-displacement curve becomes perfectly horizontal 
or vertical, several problems in solution procedure appears. 

As an example, in Figure 5.1(a), a load-displacement curve is shown. For simple UL solution algorithm 
used, a load step increment is used, and once the solution has reached point A on the curve, after the 
next load increment, however small it may be, it will jump to point C, instead of going to point B to 
capture the curve perfectly. To resolve this problem, a displacement control method can be used. 


Displacement Displacement 

(a) (b) 

Figure 5.1 : Various load - displacement curves: (a) Snap-through, (b) Snap-back 

The most simple displacement control is just incrementing the displacement field by some admissible 
amount and calculating the corresponding force to balance the system. ' 

'This process has a problem for multiple variable problems. It is next to impossible to foretell the admissible displace- 
ment fi eld the structure is going to pass through in the next time step. Assigning some specifi ed fi elds may ruin the solution 
path completely and take the converged result to somewhere else tlian the actual one. 

But for single variable problem, this works fi ne as any assumed displacement is an admissible one. 
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But this method also fails for the case shown in Figure 5.1(b), i.e a snap-back behavior of the structure. 
Then also, the solution simply jumps from point A to C ignoring the actual curve going through point 
B. 

The best solution to this problem is a mixed control technique, i.e both the displacement and the 
force fields are incremented in a very controlled and defined way to track the actual path of the load- 
displacement curve. The most elementary type of this technique is the Arc-length control method. 


5.1 Generalized load-displacement control 

This section is targeted to cover the elementary idea of various continuation methods. In general load- 
displacement control, the target is to trace the actual path of the load-displacement curve by judiciously 
controlling both the load and the displacement fields. This is done by considering displacement as 
well as load, as variable quantities in the momentum balance equation used. To solve for the extra 
unknowns, some specific relation between the displacement and the force field is used. Arc-length 
control is one of such relations. 

Here, the momentum balance is written as [9] 

r(u, A) = fint{u) - Xfext = 0 (5.1.1) 

where fext is the load pattern applied externally (A is multiplied with it to define a proportional load) 
and fint is the internal forces caused by the displacement field u. A is called the load level parameter 
and is considered as a variable along with Ui, i = 1 to n. Hence the problem becomes a (n -f- 1) 
dimensional problem to be solved with Equation 5.1.1 defining only n number of equations. 

The (n -I- equation is obtained from some relation between Uj-s and A. For arc-length control, this 
equation turns out to be [9] 

a = (Au^Au -f AA^^VL/ext) - = 0 (5.1.2) 

where A1 is the fixed length of the arc drawn from the current point on the load-displacement curve to 
get the next equilibrium point on it as shown in Figure 5.2. Au is the total increment in the displacement 
field and AA is the corresponding increment in the load level parameter. The corresponding iterative 
increments are expressed by 6u and 6X. The parameter is the a scaling parameter between force and 
displacement. For spherical arc-length method, ^ is assumed to be 1. Equations. 1.1 and Equation 5.1.2 
can be expanded by Taylor’s series about their initial equilibrium configuration 0 and using Equation 
5.1.1, as 

f (u, A) = fo -f ~ + K5u — Xfext = 0 and ^ 

a = Uo + '2Au^ Su -4- 2AA5A^'^/^t/eit = 0 
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Force 


B 



Figure 5.2; Spherical arc-length method 


Combining Equation 5. 1 .3 in a matrix form. 


K -fext / \ “ ^0 

2Au^ J I I - ^0 


(5.1.4) 


Equation 5. 1 .4 is often called the bordered equation and can be used directly to find 6u and SX. 


5.2 Spherical arc-length method 

In spherical arc-length control method, the basic idea is that the constraint surface is made to be a 
spherical hyper-surface in (n-t-l) dimensional plane around the current solution point. The predictor 
solution as well as the corrector steps for that predictor should be arc-length controlled [10]. 

Often, the load-displacement scaling parameter is taken to be one such that the constraint equation 
represents the equation of a sphere in (n-i-1) dimension as 


Au^Au + AX^fJ^Jext = 


(5.2.1) 
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This constraint equation, along with the iterative incremental momentum equation, is arranged in the 
matrix form as Equation 5.1.4 and solved for the Su-s and SX. Then the next increment step is calculated 
as 


Aui+i — Aui + 6ui 
— AAj + SXi 


(5.2.2) 


One of the problems encountered in spherical arc length method is that the constraint surface in Equa- 
tion 5.2.1 is quadratic in A A and has two possible roots for it satisfying the momentum balance as 
shown in Figure 5.3. The correct one has to be chosen to prevent doubling back of the solution. This 
is done by keeping track of the angle of the old increment and comparing the direction difference of 
it with the two possible solutions from constraint equation. The one with minimum angle with the 
previous increment can be chosen to prevent doubling back of the result [10]. This situation is shown 

AUq 



in Figure 5.3. For the case shown, the solution {(^Ui, (5Ai} is the solution which ensures no doubling 
back as it has the minimum angle with the previous increment Au direction. 

But in case of snap-back, the doubling back of the solution is a desired phenomenon. For that case, one 
possible method to find out the correct value of ^ A, is the minimum residue criteria [10]. 

To detect that whether the solution is near to the limit point or not, the current stififiiess pEirameter may 
be used [9]. This is nothing but the ratio of the current tangent stiffness to the initial tangent stiffiiess. 
It gives a rough idea about the slope of the load-displacement curve at the current solution point. 
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5.3 Linearized arc-length method 


The spherical constraint Equation 5.1.3 can be linearized by considering ao = 0, i.e the constraint 
condition to be satisfied initially [9]. 

2Au'^6u + 2AX6XfJ^J^t = 0 


or. 


<iA = - 


AvFSu 


This constraint, if arranged in proper maimer is nothing but 


(5.3.1) 


AiF AA/L SXfl, 

Au^Au SvT'Su 

i.e. the iterative incremental solution space {(ift, (iA) is orthogonal to { Au, AA}. 

Two different corrector algorithms based on this linearized arc-length method is shown in Figure 5.4. 
Figure 5.4(a) shows the Riks-Wempner scheme. Here, the corrector solution {6u, (5A} is orthogonal to 




Figure 5.4: Linearized arc-length methods: (a) Riks-Wempner method, (b) Ramm’s method 


the predictor solution { Au, AA}. And Figure 5.2(b) shows Ramm’s method in which the predictor step 
{Au, AA} is continuously updated and the corrector solution is orthogonal to the secant change [16]. 
The momentum balance along with the linearized arc-length method can be arranged in matrix form as 


^ -/ext f 

AxF AXfZtLt J 1 


- ro 

0 


(5.3.2) 
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5.4 Cylindrical arc-length method 

This method was introduced by Crisfield in 1981 [9, 16]. He modified the spherical arc-length method 
by putting 'T = 0 in Equation 5.1.2. That makes the constraint surface a cylindrical one with its axis 
along the direction of load parameter A as 


AifAu = Al^ (5.4.1) 

While controlling the arc-length keeping a specific displacement component fixed instead of the whole 
displacement field, this method reduces to the simple displacement control method and has the disad- 
vantage mentioned at the starting of this chapter. 

But this method is a very simple one and can be readily implemented in any computer program. The 
problems of line search, finding the feasible solution between etc are not issues as they are in spherical 
arc length control. Hence, this method is used to predict limit point behavior of some stmcture using 
the code. 
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Chapter 6 
Discretization 


While discretizing the equations by Finite Elements and Finite Difference, the material co-ordianate 
is considered for spatial discretization. The displacement field u(x, t), being the primary variable, is 
expressed in a semi-discrete manner by separation of time and space variables as [4, 3] 


u{x, t) = ^s{x)u^{t) 
= 


( 6 . 0 . 1 ) 


where is some interpolation function and suffices s and m denotes the spatial and material functions, 
respectively. 


6.1 Discretization in space 

In finite element discretization, a Lagrangian mesh is considered, and the corresponding approximation 
is taken by isoparametric assumption as 


X = i^{X)x%t) 
u = ^{X)u%t) 


( 6 . 1 . 1 ) 


where $(X) is the shape function matrix. 

Flence, the virtual displacement field can be discretized as 

6u = l(X)5u^ (6.1.2) 


The incremental parameters as 

Au = l(X)Au" (6.1.3) 
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and time derivatives as 


ti = ^{X)t 

ii = 

Au = $(X)AS® 

AH = l•(X)A{i® 

and the deformation gradient is then given as 

p=^ 

dX 

-- dkx).. 

dX 

6.1.1 Dynamic forces 

From Equation 6.1.1 and Equation 6.1.4, the dynamic part of Equation 3.2.10 becomes 

/ ptAu ■ 5udQ = / ptSu^AudQ, 

■JBt JBt 

= I 

JBt 


Su^'^if ptl'^l^d^yAe 

JBt 


i.e 


/ 

Jb, 


Pt 


Au ■ 5udQ = Su^'^MAu*^ 


where M = pti^^^dO. is known as the mass matrix. 


6.1.2 External forces 

From Equation 3.2. 10, the ex.temal body force part can be discretized as 


I ptAb ■ SudQ — / ptdu^Abdfl 
JBt JBt 


= / ptdu^^ ^Ab^dCt 
JBt 

= Su^'^{[ pt^^^Ab'^}dQ. 
JBt 


I.e 


I PtAb ■ 6udQ, — Su"'^ fb 

JBt 


(6.1.4) 


(6.1.5) 


( 6 . 1 . 6 ) 


(6.1.7) 
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where fh — pt^^^A.b^dil is the discretized external force vector. 
And the surface force 


AT • 6udS 


'dBt 


L 


Su^AfdS 


dBt 


/dBt 




i.e 



5udS = 5u^^fs 


( 6 . 1 . 8 ) 


where $5 is the shape function matrix for the boundary elements and is the 

discretized surface force vector. 

Hence the total external force is discretized as 


[ ptAb ■ 5udQ + [ AT- SudS = Su^^{fb + fs) = Su^'^{Af) (6.1.9) 

JBt JdSt 


where A/ = /* + /«, is the total incremental external force. 


6.1.3 Internal forces 

The incremental stress tensor in Equation 3.2.10, expressed by the Jaumann-Zaremba rate increment 
mentioned in Section 3.3, is 


A^ 



+ Wdt + 

m 


= ~Eh-V + dtW'^ 


First, the incremental part of d for stretch is considered. The corresponding virtual work is 

Ad = Ad : 5{b + W) 
ox 

= Air : SD (as W is skew and Ad is symmetric) 

= {bb + Wdt + dtW^) : sb 

That implies [10] 

[ Ad : ^da = f ibb :6b + 2dtW^ : <5£))dfl' (6.1.10) 

JBt 9x Js, 

For discrezitation, the two matrices need to be expressed in Voigt notation, or in other words, in a vector 
form such that the balance equations remain the same. 

*if a, 6 eTZsym and w eTZskew, then wa:b = tr{aw'^b) = — tr{awb) = — aifcWfcjfty = Wikakjbij = aw^ : b 
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Conventionally, the stress matrix 


^XX 

^xy 

^xz 

^yx 

^yy 

O'yz 

f^zx 

(^zy 

^zz 


is written in vector form by the Voigt rule for symmetric kinetic tensors[4] as 


{^voigt} ^yy ^zz ^yz ^xz ^xy} 


( 6 . 1 . 11 ) 


And from Viogt rale for second order symmetric kinematic tensor, the deformation tensor 


Dxx 

^xy 

Dxz 

Dyx 

^yy 

Dyz 

Dzx 

Dzy 

Dzz_ 


is arranged in vector form as 

{D^oigt} = {DxX Dyy D^z ^Dyz 2Dxz 2Z)ly}^ (6.1.12) 


such that, 


(7 . T) — ^voigt * 


voigt 


Now, the as the spin tensor is not a symmetric one, the conventional Voigt rules cannot be applied to 
transfer them into vector notation. Rather, the corresponding virtual work is expanded as follows. 


: 5D = tr 
= tr 


W^6D) 

^sb)W 


{5D)a : W 


(6.1.13) 


Now, say {5D)at = 


On 

ai2 

Ol3 


■ 0 

a P 

021 

0-22 

O 23 

,and 1V = 

— a 

0 7 


032 

033_ 


.-0 

0 

1 


(6.1.14) 
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Then 



Uii 

ai2 

ai3 


0 a 

p 


[SD][at] : = 

<^21 

^^22 

0-23 


- 

a 0 

7 



_fl31 

<^32 

®33_ 


- 

P -7 

0 _ 



1 

/ 

Oil 

0-21 

^31 


■ 0 

a 

P 

= tr 


“12 

0-22 

<3'32 


— a 

0 

7 



V 

0-13 

0-23 

033_ 


-p - 

- 7 

0 _ 


— < y { ai 2 - a 2 i ) + Pio-n ~ <2.31) + 7(023 — 032) 
Combining this result with Equation 6.1.13 and 6.1.14, 


: SD = 


(012 — 021) 
■* (fll 3 ~ Q31) 
^ (023 — 032) 



(6.1.15) 


Equation 6.1.15 gives the required Voigt vector expression for the rotational part of Equation 6.1.10. 
Now, expressing the terms by the shape function [$], deformation part of Equation 6.1.10 



kD-5bd?t= I hsLfbdQ 

JBt 


7 . 


ESu^'^B'^Bu^dn 


where E is the material tensor expressed in matrix form by Voigt notation as discussed in Section 3.6 
and B is defined as b = BAu®. Then, the internal force for deformation can be discretized as 



Eb ■ 6bdn = 


(6.1.16) 


where Kq = bb^BdO, is known as the stiffiiess matrix for the deformation b. 

And using Equation 6.1.13,6.1.14 and 6. 1 . 1 5, the rotational part of Equation 6. 1 . 1 0 is discretized as 


: 5bdn = 


f (ai2 ~ 021)'] 

(Qi3 ~ Q 31) / 

(^23 “ O 32 ) J 


iBt 


dO, 


where Bj is defined as {(ai2 — a2i), (cis — Usi), (023 — 032)}^ = Bi6u^, and B 2 as {a, /?, 7 }^ = 
J52Au®. 
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Hence, the rotational part is 



; 5Ddn = Su^^Kw^u^ 


(6.1.17) 


where Kw = Bj B^dQ., is the stiffness matrix for the rotational parts, or the geometric stiffness 
matrix of the structure. 

So, the total virtual work done by the internal forces is 



ox 


(6.1.18) 


where K = + Kw is the total stiffness matrix of the system. 


6.1.4 Expressions of different matrices used for discretization 

Say, the primary unknowns {wi, U2, us} at any generic point is interpolated as 


Ui = + ... + NnnU\ "" 

U2 = Niuf + Niuf + ... + NnnUl "" 
and Us = Niu^ + IS/ 21^^ + ... + Nnnu\ 


Or, 


M 

>1 

0 

0 

N 2 

0 

0 ... 

- Nnn 

0 

0 ' 

L U 

0 

iVi 

0 

0 

N 2 

0 ... 

... 0 

Nnn 

0 

U 

_0 

0 

iVi 

0 

0 

to 

... 0 

0 

Nnn_ 


u' 


u: 


el 


U 


yC nn 

/i 7%Tl 

U2 


Hence the shape function matrix is given by 


[$] 


A^i 0 0 N 2 0 0 

0 ATi 0 0 iVa 0 

0 0 iVi 0 0 N2 


Nnn 0 0 

o' Nnn 0 

0 0 Nnn 


(6.1.19) 
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The rate of deformation tensor (in Voigt notation) can be expressed as 




IBKAu} 


i.e. 

$11 $12)Xl — ^1 nmxl 

^21-,x2 $22ii2 $2nn)i2 

$31)i3 $32,13 ■••• $3nn,x3 

$21 5X3 +$31 5x2 $22,x3 +$32 5x2 — ^2 nnixZ +$3 Tin 5x2 
^115x3+^315x1 5^125x3+^325x1 •••* nn,x3 +$3 nniil 

$ll,i2 +$21,xl $12,i2 +$22,il $1 nn,x2 +*^*2 nn,il 

Similarly, the other strain displacement matrix, for rotation, can be written as 

-^1,12 “■■^lixl 0 A^2 ,x2 ~-^2,x 1 0 ••• •^nn,x2 ~'-^nn,xl 0 

[• 82 ] = - Ni,x 3 0 — A^ljil N2,x3 0 —N2,x1 ••• Nnn)x3 0 ~-^nn,xl 

L 0 N, 5X3 -ATi ,i2 0 ^2ixZ ^2 1x2 0 ^nnixS -^nn)x2j 

( 6 . 1 . 21 ) 

And now combining the Equations 6. 1 .6, 6. 1 .9 and 6. 1 . 1 8, the incremental virtual work statement can 
be arranged as 

T Su^'^kAu^ = Su^'^Af 

As Su is an arbitrary virtual displacement, the discretized momentum balance is 

[M]{Au^} + [K]{Au^} = {Af} (6.1.22) 



( 6 . 1 . 20 ) 


6.1.5 Discretized heat transfer equations 

The temperature field is assumed to be interpolated as 

9 = Niei + N2e^^ + ... + NnJnn 

or, e = {$i}"’{6^} 
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Then 


se = 

M = {^>:}^{A^^} and 
A(9,, = {<I>i,J'"{A0^} 


Now, the incremental energy balance equation for heat transfer (Equation 3.5.5) can be discretized 
similar to its elastic counterpart as 


'Bt 


zp{69^y{^^}{^,}^{Ar}dn+ I {ki{6eY{Ki}{^i,y{^o^})dn 

Jst 

= [ + {■).} 
JBt 


Or, 


{■59'}"' (la»l{A^} + [ir,J{A«'}) = {59*}''{A,,} 

where |C.J = c/){4>iH'*i}’'<ia and [Kai = Js, (ii{'59'}"'{4i,i}) dfi. i = 1, 2, 3 

and {Ag} = + {gs}, {gj} being the heat flux through the boundary. 

Again, as {S6^} can be chosen arbitrarily, the discretized heat flux balance can be written as 

[Cth]{^9^} + \KtH]{^9^] = {Ag} (6.1.24) 

6.2 Discretization of the contact constraint 

This work includes just the sticking contact of the body against a rigid wall. For that case, the penetra- 
tion function for the Lagrangian multiplier formulation in Section 4.2 becomes 


9 — '^contact 


where Uccmtact are the displacement fields at the contacting nodes. These contacting nodes are a subset 
of all the nodes and the penetration function can be written as' 

{g} = [C]K} 

if Uj is the contacting d.o.f (6.2.1) 

otherwise 


where Cij = 
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And also, 

{<9} = |C){fa'} 

and {Ag} = [C']{Aw®} 

Then Equation 4.2.5 can be discretzed as 

{6uY[M]{Au^} + {6u^nK]{Au^} - {5u^Y{AU} 

- ({<5u^}^[C]^{AA} + {(5Af [C]{Au*}) = 0 

Or, 

([M]{Au^} + [K]{Au^} - {A/,,*} - [a]''{AA}) = 0 

and {<^A}^[C]{Au®} = 0 

where A is the contact forces at the contact nodes. 

As both the virtual fields {Su^} and {A} can be chosen arbitrary, the momentum balance turns out to be 

+ [K\{/\v.‘) - |C1’'{AA} = { A/„,} and 

^ (o.z.z) 

[<:7]{Au=}=0 

Equation 6.2.2 can be combined to be arranged in form as 

\m] [0]1 r {Au«}l \[K] -[C]^] ({Au^}\ ^ /{A/e^jl 

_ [0] [0]J \ {AA} / |_[C] [0] J \ {AA} j \ {0} j (6.2.3) 

or [Mcontact]{Au'} + [Kcontact] { Au^ } = {AE} 

where {An'} = {{Au^}, {AA}} and (A/'} = {{A/ext}, {0}}. 

6.3 Discretization in time 

Conventionally, two types of methods are used for time integration scheme of Equation 6. 1 .22. Namely, 
explicit and implicit time integration scheme. 

Explicit schemes are based on the idea that equilibrium vales are known at time t, and they are used to 
obtain the solution at time t+ At [4, 1 7]. As the name suggests, this scheme is a rather direct method to 
integrate the time dynamic equation and no prior transformation is done before integration [3]. Though 
the explicit integration schemes like Central difference scheme, are very popular and simple to apply, 
they have one disadvantage that they are conditionally stable [17]. 

The implicit schemes of numerical integration, on the other hand, are unconditionally stable. They 
depend on the values of the variable and its derivatives in both time steps t and t + At. One of the most 
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popular methods in implicit schemes is the Newmark-(3 method. This method is often interpreted as 
the extension of linear acceleration method [3]. The assumptions in this method are 


= Ut + UtAt + 


(- - a)ut + aut+At 


Ar 


(6.3.1) 


and 

Ut+At =Ut+ [(1 - S)ut + ^Ut+Af] At (6.3.2) 

where a and 5 are Newmark constants whose values dictates the dependency between the values and 
slopes at time t and t + At. In this work, the constant acceleration method, or the trapezoidal rule is 



Figure 6 . 1 ; Constant average acceleration Newmark scheme 

used for which ci = | and ^ = 5 . This scheme is shown in Figure 6.1. This method is unconditionally 
stable and was originally proposed by Newmark [3]. 


6.3.1 Modification of discretized equations for Newmark integration scheme 

From Equation 6.3.1 , 


Au 


UtAt + 


(2 ~ + OUj+At 


At^ 


or. 


And from Equation 6.3.2, 


1 1 . 1 - 

Au = — 7 -x Au — Ut — — Ut 

OiAt^ aAt 2a 


(6.3.3) 


Au = [ut + (5Au] At 


Replacing Au from Equation 6.3.3, 

Au = 


S S ^ 

Ut + Aw - -T-ut - —Ut 


aAi^ 


aAt 2a 


At 


(6.3.4) 
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Using Equation 6.3.3 and 6.3.4, Equation 6 . 1.22 can be modified as 


[K]{Au®} = {AF} 

where, [if] = [if] + ( 6 . 3 , 5 ) 

and {A/} = {A/} + + i[M]{«,} 

Similarly, the discretized heat transfer Equation 6.1.24 can be modified using Newmark relation Equa- 
tion 6.3.4 as 


[Kth]{A0"} = {Aqth} 
where, = [KtH] + 

and {/\qth} = {Ag} -f ~[Cth]{dt} - (l 

a \ 

6.3.2 Stability of Newmark-/5 scheme 

As discussed in Cook [13], the accuracy, stability and damping is also a matter of importance while 
using Newmark-/3. The scheme is said to be unconditionally stable for 

2a = 5=l- 
2 

But as observed from the computational results, often there exists an oscillation in the solution with 
a = 0.25 and S = 0.5. This oscillation can be reduced using higher values of those constants. But 
according to Cook, the process reduces the accuracy of the solution. For (i > 0.5, the stability is 
conditional and the condition on the corresponding time step length is [13] 

^ ^critTmin Q; < ^ ^ 

W-max 2II 2 

where Wmax and Tmin are the highest natural frequency of the structure and corresponding time period. 
And Clcrit is defined as 


(6.3.6) 


2q; 


^t[Cth]{0t} 


" 5/2 -a 

where ^ is the damping ratio. When 6 > 0.5, allowable time step length is increased damping. 

For (5 > 0.5, Newmark method displays algorithmic damping, and according to Cook, the accuracy of 
the result is reduced. In that case, the right choice of a is a = \ (<^ + |)^. 
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Chapter 7 

Results and discussion 


For checking the theories deduced, a computer program is developed using FORTRAN 90. 
Tri-quadratic Lagrangian shape functions are used for spatial finite element discretization and constant 
acceleration Newmark-/? or trapezoidal time integration rule is used for discretization in time. Newton- 
Raphson scheme is used for the predictor-corrector steps used. 

As a first step, some standard and benchmark static problems are solved by UL method in a quasi-static 
maimer (refer Figure 2.7). And for thermo-elastic problems, heat transfer equations are assumed to be 
decoupled from the elastic equations. First, heat transfer problem is solved independently and then the 
resulting stresses are fed into the elastic part to solve for the deformation. Then, a few of the continu- 
ation techniques are coded to capture some standard unstable phenomena. The contact constraints are 
also applied and their validity is checked. 

But for the dynamic formulation, Newmark-/? time integration scheme is used which is unconditionally 
stable for specific values of the Newmark constants. The condition for the process being uncondition- 
ally stable is 2tt > 5 > | [13]. And for conditional stability, the time step size should be less than 
some critical size [13]. 

For the inherent instability and slow convergence rate of Newmark-/? scheme, the continuation meth- 
ods have not been extended to the dynamic cases as they themselves are quite complex regarding the 
application of Newton-Raphson method near the limit points. 

The results have been compared with strength of material solutions and ANSYS results using 20 noded 
brick elements. The material constants of the structure is taken as shown in Table 7. 1 


Young’s modulus 

2.1 

Poisson’s ratio 

0.3333 

Density 

7.8 

Thermal conductivity 

0.2 

Thermal expansion coefficient 

1.0 e-® 


Table 7.1; Material constants used 
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7.1 Static analysis 


The model used for the static analyses is mainly a uniform rectangular cross section beam / bar as 
shown in Figure 7. 1 . For the bending and axial tension, the cross section of the beam is taken as 1 unit 
X 1 unit and the length is taken as 5 units. For the buckling analysis, the length is taken to be 20 units. 





Figure 7.T. Model used for static analysis 


7.1.1 Bending of a beam 

A distributed flexure load is applied at the free end of the beam models shown. The results, compared 
to strength of materials(S.O.M) and ANSYS are as follows: 


End load 

End deflection 
from S.O.M 

End deflection 
from ANSYS 

End deflection 
from Program 

10 units 

2.3810 e-2 

2.4027 e-2 

2.4011 e-2 

20 units 

4.7619 e-‘^ 

4.8054 e-'-^ 

4.7880 e-^ 


Table 7.2: Load vs deflection for bending of a beam for static analysis 


From Table 7.2, it can be observed that both the solutions obtained by ANSYS and the program are to 
2 some extent softer than the one obtained by strength of materials Euler Bernoulli beam theory. This 
difference is expected for the constraints imposed on bending by Euler’s theorem. 

But though the effects of rotation is considered in the program (through Jaumann-Zaremba objective 
stress rate) the ANSYS results happen to be a little bit softer than the program output. Apparently, this 
is because of the comer singularities of the beam eind the difference in the types of elements used. The 
code uses 27 noded tri-quadratic Lagrangian element while in ANSYS, the closest element available is 
the 20 noded serendipity brick elements. 

The load-displacement curve for this results is shown in Figure 7.2. 
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Figure 7.2: Load-deflection curve for bending load 

7.1.2 Axial end load applied on a bar 

The same model is used as before. A uniformly distributed tensile load is applied at the free end. The 
corresponding deflections are found as: 


End load 

End deflection 
from S.O.M 

End deflection 
from ANSYS 

End deflection 
from Program 

50 units 

1.1905 e-^ 

1.2040 e-3 

1.2699 e-^ 

100 units 

2.3810 e-^ 

2.4080 

2.4930 e-“ 


Table 7.3: Load vs deflection for axial stretch of a bar for static analysis 


Table 7.3 shows a clear difference between the strength of materials. O.M), ANSYS and program 
solutions for the deflection of the free end under axial tension. The ANSYS solution is a bit softer than 
the s.o.m solution, and the solution obtained from the program using 10 updated Lagrangian steps is 
much more softer. This kind of result is expected as updated Lagrangian scheme was used and effect 
of rotation was taken into account. The plot of these results is shown in Figure 7.3. 






Figure 7.3: Load-deflection curve for axial load 


7.1.3 Heat transfer analysis 

The thermal balance equation is considered to be decoupled from the elastic one in this analysis. The 
model used in the previous sections, is taken and a pure thermal load (no force or moment) is applied. 
The temperature response at the free end of the bar under a uniformly distributed end heat flux is 
shown in Table 7.4. Here also, the updated Lagrangian steps gives a high numerical value of the free 
end temperature. 


End heat flux 

Free end temperature 
from ANSYS 

Free end temperature 
from Program 

10 units 

250.0000 K 

255.6780 K 

20 units 

500.0000 K 

511.3561 K 


Table 7.4: Temperature at the end of a bar with uniform heat flux at the free end 


7.2 Application of continuation techniques (static analysis) 

The continuation methods, mentioned in Chapter 5, are applied on some standardized models and are 
compared with the compared with conventional results. The most simple continuation algorithm, the 
cylindrical arc-length algorithm [9] is used in the program to capture the structural instabilities. 
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7.2.1 Buckling of a beam 

The same model mentioned above is used with a compressive end load. The length of the bar is taken 
to be 20 units, or rather L/D ration is taken to be 20. As the non-linear balance equations has been 



Figure 7.4: Model used for buckling analysis 

linearized about the equilibrium configuration, a very small perturbation flexure load should be applied 
at the end to initiate buckling as shown in Figure 7.4 

In Figure 7.5 the result from the cylindrical arc length method is plotted against the critical buckling 
load obtained from Eulers formula for long columns. 



Figure 7.5: Load-displacement curve for buckling of a slender bar 


L/D ratio 

Critical Euler load 

Critical load from program 

20 

1727.1808 1717.4666 (at 0.1 unit end deflection) 


Table 7.5: ANSYS and program solution for buckling load 
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In the load deflection curve shown in Figure 7.5, the critical load obtained by the computational method 
can be taken as the one where the curve becomes horizontal just before drooping down (at displacement 
0.1 approximately). 

7.2.2 Snap through phenomenon in a shallow arch 

Snap through phenomenon is a unstable behavior of a structure under loading. It is a situation where 
load parameter decreases with an increase in displacement. This phenomenon is shown in Figure 5. 1 (a). 
A shallow arch model shown in Figure 7.6, is tJiken for the analysis. A central load is applied at the 



Figure 7.6; Model used for analysis of snap-through behavior 

middle of the arch and corresponding load-deflection curve is shown in Figure 7.7. The cylindrical 
arc-length method is able to capture the snap-through behavior, but the conventional FEM solution is 
not. Rather, shown by dashed line, it almost jumps through the drooping portion of the curve as shown 



Figure 7.7: Comparison of cylindrical arc-length method with conventional solution method for snap- 
through phenomenon in a shallow arch 

schematically in Figure 5.1(a). But here, the curve resulting from conventional FE solution lies some- 
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ANSYS results are generated using displacement boundary constraints at those nodes are close to the 
values obtained by the program. 

The deformed shape along with the undeformed edges of the model used is shown in Figure 7.10. 



Figure 7.10: Deformed shape and undefonjied edges from static contact analysis 
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7.4 Dynamic analysis 


For dynamic analysis, the differential equation in time is discretized using Nemark-/3 algorithm. Ini- 
tially the trapezoidal version of Newmark algorithm is used. 

7.4.1 Bending analysis - stability of Newmark method 

Newmark-/? method with a = 0.25 and S = 0.5 is applied for the beam model under flexxne loading as 
shown in Section 7.1.1 with zero initial conditions. The results obtained are shown in Table 7.7. 


End load 

ANSYS result 

Program with 5 
UL steps 

Program with 10 
UL steps 

Program with 15 
UL steps 

1 0 units 

2.4029 e-2 

2.3958 e-2 

2.6067 e-2 

2.0019 e-2 


Table 7.7: End deflection for a beam under flexure load - dynamic analysis 


These results, when plotted as load-displacement curves (Figure 7.11), an oscillation is observed in 
the solution. The oscillation is more for UL steps with smaller time step lengths. This oscillation is 



Figure 7.1 1: Increasing oscillations in Newmark-/? technique with decreasing time step length using 
Newmark constants a = 0.25 and S = 0.5 (trapezoidal rule) 

seemingly because of the low values of the Newmark constants {a = 0.25, S = 0.5) taken. The values 
of the Newmark constants, though mentioned as [13] 


2a > <5 > - 

jL 
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for unconditional stability, there remains some restrictions on the step lengths as seen in Figure 7. ! 1. 
As suggested by Mahato, the numerical values of the Newmark constants are increased (to a = 1.56 
for S = 2.0, from the condition mentioned in Section 6.3) to get a stabilized load-displacement curve 
as shown in Figure 7.12. 

This increment of the numerical values of the constants introduce algorithmic damping [13] to damp 



Figure 7.12: Stabilized Newmark-/(3 method with algorithmic damping (a = 1.56 and 5 = 2.0) 

out the oscillations shown in Figure 7.11 and then, for conditional stability, some restriction on the time 
step lengths has to be imposed. This condition is given in Section 6.3. 

7.4.2 Application of contact conditions 

To damp out the oscillations in the solution, the values of the Newmark constants are kept as q = 1.56 
and S = 2.0 for dynamic analysis. With these values, the contact conditions are tested for the dynamic 



case for another model. The model is a beam, supported on hinges at two ends and loaded at the middle 
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as shown in Figure 7. 13. The hinged joints are simulated by rigid contact for the corresponding nodes. 
A distributed force is applied on the middle line vertically downwards and corresponding displacement 
is computed. 

The results are given in Table 7.8. 


Load 

ANSYS result 

Program solution 

10 units 

4.3434 e-2 

5.0806 e-2 

20 units 

0.08686 

0.12541 


Table 7.8: Mid point deflection of a hinged hinged beam - dynamic contact analysis 

The results obtained from the code in this case, unlike the static contact analysis, is more flexible 
as compared to ANSYS. This is probably because of the algorithmic damping used by the increased 
numerical values of the Newmark constants. The deformed shape of the model is shown in Figure 7. 14 



Figure 7.14: Deformed shape and undeformed edges from dynamic contact analysis 


7.5 Thermoelastic deformation (dynamic analysis) 

The model used for static analysis with axial force, is used for this case. One side of the bar is fixed 
and some force as well as heat flux is applied on the free end. 

The results obtained from the program and ANSYS is listed in Table 7.9 


End load 

End heat 
flux 

End deflection 
from ANSYS 

End temp from 
ANSYS 

End deflection 
from program 

End temp from 
program 

100 

50 

4.958 e-^ 

2545 K 

4.8161 e-^ 

2517.43 K 

200 

50 

9.915 e-^ 

2545 K 

9.6398 

2519.36 K 


Table 7.9: Dynamic thermoelastic analysis 
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As the thermal analysis in ANSYS was carried out independently, the temperature for both the cases 
are same. The elastic response obtained from the simulation yet remains a little stiffer than ANSYS. 
For almost all the analyses simulated by the program is a little stiffer than the solution obtained through 
ANSYS. The existance of comer singularities in the models used may be responsible for that and a 
graded or adaptive mesh may be proposed for the accurate analysis of the stmctures. 
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Chapter 8 


Conclusions and futnre scope 


8.1 Conclusions regarding the present work 

The target of the study is mainly to develop a formulation for the much discussed updated Lagrangian 
scheme. The basic definition of the method is used along with transformation laws from continuum 
mechanics along with proper approximation to get a honest incremental form of the UL method. In the 
literatures available, there are some small deformation assumptions plugged into the UL formulation. 
The general formulation without those approximations may be really useful for large strain-large rota- 
tion problems. Also the step size used in UL method becomes a parameter that one can control. With 
fewer assumptions, it is possible to increase the step sizes reducing the computational cost drastically. 
Thermodynamics is also included in the code to widen its possible scope of application. 

In deducing the incremental form, objectivity of physical kinematic quantities are taken care of to make 
the analysis able to handle problems with large rigid rotation. Oldroyed rate, a general objective mea- 
sure of stress rate is discussed, and the Jaumann-Zaremba rate form, a derivative of it and one of most 
popular stress rates, has been used in the computer program. 

The contact formulation has been developed for both sticking and slipping contacts. The code is de- 
veloped for only sticking contact with a rigid wall, but the Lagrangian multiplier technique is used in 
a flexible manner such that the frictional contact models and the two body contact problems can be 
modeled afterwards. 

Some continuation techniques have been reviewed with the target to make the solution scheme more 
efficient. One of the continuation methods, the cylindrical arc-length method, has been plugged into 
the code to test the formulations used and a few global instable phenomenon have been captured. 

After all, for the dynamic analysis, a problem regarding the oscillating nature of the solution fi’om 
Newmark-/? time integration scheme is faced. Some literatures have been referred regarding this issue, 
and following those, the oscillations had been damped out using algorithmic damping. 

Finally, some standard models are used to validate the code. The problems, though very simple in 
nature, shows the scope of application of the code and its effectiveness. The numerical values obtained 
from it are compared with strength of material and/or ANSYS solution. The results obtained, expected 
to be a bit flexible than ANSYS, is not exactly satisfactory for all the cases. But they remain very close 
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to the bench-mark solution and ANSYS results. 


8.2 Scope of further work 

The idea behind developing the incremental updated Lagrangian method from the basics was mainly to 
get a general expression for the same without any pre-assumed points. This form seemingly can enable 
one to follow almost the actual path in solving a problem and the time steps of the solution procedure 
can be reduced. But as the current analysis is a linear one, the rejection of the higher order terms has 
put a limit on the number of UL steps used. A non-linear analysis and corresponding simulation may 
be carried out to focus on the reduction of the computational efforts. 

The assumption discussed in Section 3.2, that the virtual displacement at time t and t + At are same, is 
not always valid, specially for large displacement contact problems. The virtual displacement Su itself, 
may be considered as a function of the space and time which will lead to a more general incremental 
weak form. 

Also the application of the contact constraints has been a very restricted field in the present study. The 
present code may be extended to handle frictional contact and multi-body contact problems. 

Among the continuation methods discussed, only the cylindrical arc length method has been success- 
fully incorporated in to the computer program. But as mentioned in Section 5.4, this method is a cmde 
one for predicting the correct load-displacement curve followed. The spherical as well as the linearized 
arc length methods may be included in the code and a comparative study may be carried out among 
them. 

Another important issue regarding the set-backs of the simulation is the Newmark time integrator 
scheme used. Originally, trapezoidal scheme proposed by Newmark was used. Afterwards, the New- 
mark constants used have been modified (their numerical values increased) to stop the oscillations 
shown by the original scheme. This procedure seemingly induces algorithmic damping in the system. 
This may make the results to deviate from their actual values even though the spatial formulations have 
been carried out with precision. Some better time integration scheme should replace the Newmark one 
in the analysis. 

Lastly, the possible issues, like comer singularities, causing the computational result to be a stiffer one 
though large strain and deformation has been considered, may be taken care of by a good meshing and 
solver algorithm. This will make the code a more dependable and robust one. 
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